******************************
* Stata preamble
******************************
clear all
macro drop _all
capture log close
set linesize 255
set mat 1500
set more off,perm

				
cd /*set to Directory with BW_RCFS_Replication_Dataset.dta*/
capture log close
log using log_02_RCFS_Specification_Curve.txt, replace text
use BW_RCFS_Replication_Dataset.dta , clear

******************************************
* Perform Estimation and Save Results
******************************************
gen z_Mcaid_x_BigBoost = z_Mcaid_share * Mcaid_Boost

file open myfile using "SpecCurve_tab_sep.txt", write replace
file write myfile  "Model" _tab "TargetVariable_1" ///
                  _tab  "beta_1"   _tab  "standarderror_1" ///
				  _tab "TargetVariable_2" ///
                  _tab  "beta_2"   _tab  "standarderror_2" ///
				  _tab "TargetVariable_3" ///
                  _tab  "beta_3"   _tab  "standarderror_3" ///
                  _tab  "loglikelihood"   _tab  "SizeMeasure"  ///
                  _tab  "FullCommand"  _tab "Timestamp" _n

eststo clear	
local SIZE z_Beds z_Beds_2 z_Beds_3 z_AvgResidents z_AvgResidents_2 z_AvgResidents_3 z_admissions_hundreds
qui eststo: nbreg covid_cases z_logCash z_Mcare_share z_Mcaid_share z_Mcaid_x_BigBoost	`SIZE'  z_FacilityAge  z_ZipBlackShare  z_LPNhrsPerResDay z_RNhrsPerResDay z_TotalMargin z_Leverage ForProfit i.HealthInspRating i.CntyFE  ,  cluster( CntyFE )
	file write myfile	%9s "`e(cmd)'" ///
						_tab %30s 	"z_Mcare_share" ///
						_tab %12.4f	(_b[z_Mcare_share]) ///
						_tab %12.4f	(_se[z_Mcare_share]) ///
						_tab %30s 	"z_logCash" ///
						_tab %12.4f	(_b[z_logCash]) ///
						_tab %12.4f	(_se[z_logCash]) ///
						_tab %30s 	"z_Mcaid_x_BigBoost" ///
						_tab %12.4f	(_b[z_Mcaid_x_BigBoost]) ///
						_tab %12.4f	(_se[z_Mcaid_x_BigBoost]) ///
						_tab %12.4f	(e(ll)) ///
						_tab %100s	"`SIZE'" ///
						_tab %500s	"`e(cmdline)'"	///
						_tab %100s "$S_DATE $S_TIME" _n					

qui foreach SIZE in "z_Beds  Occupancy" ///
				"z_Beds z_Beds_2  Occupancy" ///
				"z_Beds z_Beds_2 z_Beds_3  Occupancy" ///
				"z_Beds z_Beds_2 z_Beds_3 Occupancy z_admissions_hundreds" ///
				"z_AvgResidents z_AvgResidents_2" ///
				"z_AvgResidents z_AvgResidents_2 z_AvgResidents_3"  ///
				"z_AvgResidents z_AvgResidents_2 z_AvgResidents_3 z_admissions_hundreds"  ///
				"z_Beds  z_AvgResidents" ///
				"z_Beds  z_AvgResidents Occupancy"  ///
				"z_Beds z_Beds_2  z_AvgResidents z_AvgResidents_2 " ///
				"z_Beds z_Beds_2  z_AvgResidents z_AvgResidents_2 z_admissions_hundreds" ///
				"z_Beds z_Beds_2  z_AvgResidents z_AvgResidents_2 z_admissions_hundreds Occupancy" ///
				"i.Beds20bins" ///
				"i.Beds20bins z_admissions_hundreds Occupancy" ///
				"i.AvgResidents20bins" ///
				"i.AvgResidents20bins z_admissions_hundreds Occupancy" ///
				"i.Beds20bins i.AvgResidents20bins" ///
				"i.Beds20bins i.AvgResidents20bins  z_admissions_hundreds Occupancy" ///
				"z_Beds z_Beds_2 z_Beds_3 z_AvgResidents z_AvgResidents_2 z_AvgResidents_3" {
	eststo: nbreg covid_cases z_logCash z_Mcare_share z_Mcaid_share z_Mcaid_x_BigBoost	`SIZE' z_FacilityAge  z_ZipBlackShare  z_LPNhrsPerResDay z_RNhrsPerResDay z_TotalMargin z_Leverage ForProfit i.HealthInspRating i.CntyFE,  cluster(CntyFE)
		file write myfile	%9s "`e(cmd)'" ///
							_tab %30s 	"z_Mcare_share" ///
							_tab %12.4f	(_b[z_Mcare_share]) ///
							_tab %12.4f	(_se[z_Mcare_share]) ///
							_tab %30s 	"z_logCash" ///
							_tab %12.4f	(_b[z_logCash]) ///
							_tab %12.4f	(_se[z_logCash]) ///
							_tab %30s 	"z_Mcaid_x_BigBoost" ///
							_tab %12.4f	(_b[z_Mcaid_x_BigBoost]) ///
							_tab %12.4f	(_se[z_Mcaid_x_BigBoost]) ///
							_tab %12.4f	(e(ll)) ///
							_tab %100s	"`SIZE'" ///
							_tab %500s	"`e(cmdline)'"	///
							_tab %100s "$S_DATE $S_TIME" _n					
noisily disp "$S_DATE $S_TIME"

}

file close myfile

******************************************
* Figure IA.2: Specification Curves: Controlling for Size
******************************************
	******************************
	* Medicare Share
	******************************
	clear
	import delimited "SpecCurve_tab_sep.txt"
		gen main = (_n==1)

		gen irr_beta_1 = exp(beta_1)
		gen irr_standarderror_1 =  exp(beta_1) * standarderror_1
		gsort beta_1	
		gen n = _n	
			gen lb90 = irr_beta_1 - 1.645 * irr_standarderror_1
			gen ub90 = irr_beta_1 + 1.645 * irr_standarderror_1
			
			su lb90,d
			local mm `r(min)'
			local m `r(max)'
			local N `r(N)'
		twoway rcap ub90 lb90 n , lcolor(gs12) lwidth(medium) ///
				ytitle("coefficient on" "zMedicareShare", width(25) orientation(hor)) ysc(range(.8 1.2)) ylabel(.8(.05)1.2) yline(1) xlabel(1(1)`N',labsize(small)) xtitle("") ///
				plotregion(margin(zero)) graphregion(margin(l+3 r+5)) ysize(2) xsize(6)  ///
			|| scatter irr_beta_1 n, msymbol(O) mcolor(gs12) ///
			|| rcap ub90 lb90 n if main==1, lcolor(blue) lwidth(medium) ///
			|| scatter irr_beta_1 n if main==1, msymbol(O) mcolor(blue) ///
			title() legend(off) name(Mcare90,replace) 
		
		*****************************************
		* Get the unique values of covariates
		*****************************************
		preserve
		capture drop fc2split 
		gen fc2split = stritrim(strtrim(fullcommand))
		split fc2split , parse(" ")
		drop fc2split
		reshape long fc2split, i(n) j(j)
			keep fc2split
			drop if inlist(fc2split," " , "," , ")" , "cluster(" , "cluster(CntyFE)" , "CntyFE" , "nbreg" , "covid_cases" , "z_Mcare_share")
			bysort fc2split: keep if _n==1
			sort fc2split
			replace fc2split = subinstr(fc2split, ".", "_", .)
			levelsof fc2split ,clean
			global LEV `r(levels)'
		restore
	   foreach l in $LEV {
			noisily disp "`l'"
			  capture gen `l'  = strpos(fullcommand,"`l' ") >0 
	   }
	   foreach fe in CntyFE Beds20bins AvgResidents20bins {
			capture replace i_`fe' = 1 if  strpos(fullcommand,"`fe'")>0	
	   }
	   rename z_Beds Beds
	   rename z_Beds_2 Beds2
	   rename z_Beds_3 Beds3
	   rename i_Beds20bins Beds20bins
	   rename z_AvgResidents AvgResidents
	   rename z_AvgResidents_2 AvgResidents2
	   rename z_AvgResidents_3 AvgResidents3
	   rename i_AvgResidents20bins AvgResidents20bins
	   rename z_admissions_hundreds Admissions
	   rename Occupancy Occupancy
	   rename i_CntyFE OtherControls
		global V_LIST Beds Beds2 Beds3 Beds20bins AvgResidents AvgResidents2 AvgResidents3 AvgResidents20bins Admissions Occupancy OtherControls
		keep $V_LIST
	   *Save as matrix
	   mkmat  $V_LIST , matrix(VarPresent)
	   matrix Vt = VarPresent'
	   matrix list Vt
		heatplot Vt, colors(mono,reverse) legend(off) xlabel(,labsize(small)) ytitle("Included Variables",width(0)) ylabel(,  labsize(small)) ysize(2) xsize(6) name(vars, replace) plotregion(margin(zero)) graphregion(margin(t-5 b-12))

		graph combine Mcare90 vars, cols(1) ysize(3) xsize(6) name(Mcare,replace)
			graph drop Mcare90 vars


	*********************************
	* Cash 
	*********************************
	clear
	import delimited "SpecCurve_tab_sep.txt"
	gen main = (_n==1)

	gen irr_beta_2 = exp(beta_2)
	gen irr_standarderror_2 =  exp(beta_2) * standarderror_2
		gsort beta_2	
		gen n = _n	
			gen lb90 = irr_beta_2 - 1.645 * irr_standarderror_2
			gen ub90 = irr_beta_2 + 1.645 * irr_standarderror_2
			
		su lb90,d
		local m `r(min)'
		local N `r(N)'
	twoway rcap ub90 lb90 n , lcolor(gs12) lwidth(medium) ///
			ytitle("coefficient on" "zlog(Cash)", width(23) orientation(hor)) ysc(range(.8 1.2)) ylabel(.8(.05)1.2) yline(1) xlabel(1(1)`N',labsize(small)) xtitle("") ///
			plotregion(margin(zero)) graphregion(margin(l+3 r+5)) ysize(2) xsize(6)  ///
			|| scatter irr_beta_2 n, msymbol(O) mcolor(gs12) ///
			|| rcap ub90 lb90 n if main==1, lcolor(blue) lwidth(medium) ///
			|| scatter irr_beta_2 n if main==1, msymbol(O) mcolor(blue) ///
			title() legend(off) name(z_logCash90,replace)
		noisily disp "z_logCash"	
		gtoplevelsof sizemeasure if ub90<0
			
		*****************************************
		* Get the unique values of covariates
		*****************************************
		preserve
		capture drop fc2split 
		gen fc2split = stritrim(strtrim(fullcommand))
		split fc2split , parse(" ")
		drop fc2split
		reshape long fc2split, i(n) j(j)
			keep fc2split
			drop if inlist(fc2split," " , "," , ")" , "cluster(" , "cluster(CntyFE)" , "CntyFE" , "nbreg" , "covid_cases" , "z_logCash")
			bysort fc2split: keep if _n==1
			sort fc2split
			replace fc2split = subinstr(fc2split, ".", "_", .)
			levelsof fc2split ,clean
			global LEV `r(levels)'
		restore
	   foreach l in $LEV {
			noisily disp "`l'"
			  capture gen `l'  = strpos(fullcommand,"`l' ") >0 
	   }
	   foreach fe in CntyFE Beds20bins AvgResidents20bins {
			capture replace i_`fe' = 1 if  strpos(fullcommand,"`fe'")>0	
	   }
	   rename z_Beds Beds
	   rename z_Beds_2 Beds2
	   rename z_Beds_3 Beds3
	   rename i_Beds20bins Beds20bins
	   rename z_AvgResidents AvgResidents
	   rename z_AvgResidents_2 AvgResidents2
	   rename z_AvgResidents_3 AvgResidents3
	   rename i_AvgResidents20bins AvgResidents20bins
	   rename z_admissions_hundreds Admissions
	   rename Occupancy Occupancy
	   rename i_CntyFE OtherControls
		global V_LIST Beds Beds2 Beds3 Beds20bins AvgResidents AvgResidents2 AvgResidents3 AvgResidents20bins Admissions Occupancy OtherControls
		keep $V_LIST
	   *Save as matrix
	   mkmat  $V_LIST , matrix(VarPresent)
		   matrix Vt = VarPresent'
		   matrix list Vt
		heatplot Vt, colors(mono,reverse) legend(off) xlabel(,labsize(small)) ytitle("Included Variables",width(0)) ylabel(,  labsize(small)) ysize(2) xsize(6) name(vars, replace) plotregion(margin(zero)) graphregion(margin(t-5 b-12))

			graph combine z_logCash90 vars, cols(1) ysize(3) xsize(6) name(z_logCash,replace)
				graph drop z_logCash90 vars

	****************
	* Mcaid Boost
	*****************
	clear
	import delimited "SpecCurve_tab_sep.txt"
	gen main = (_n==1)

	gen irr_beta_3 = exp(beta_3)
	gen irr_standarderror_3 =  exp(beta_3) * standarderror_3
		gsort beta_3	
		gen n = _n	
			gen lb90 = irr_beta_3 - 1.645 * irr_standarderror_3
			gen ub90 = irr_beta_3 + 1.645 * irr_standarderror_3
			
		su lb90,d
		local m `r(min)'
		local N `r(N)'
	twoway rcap ub90 lb90 n , lcolor(gs12) lwidth(medium) ///
			ytitle("coefficient on" "zMcaid*Boost", width(23) orientation(hor)) ysc(range(.8 1.2)) ylabel(.8(.05)1.2) yline(1) xlabel(1(1)`N',labsize(small)) xtitle("") ///
			plotregion(margin(zero)) graphregion(margin(l+3 r+5)) ysize(2) xsize(6)  ///
			|| scatter irr_beta_3 n, msymbol(O) mcolor(gs12) ///
			|| rcap ub90 lb90 n if main==1, lcolor(blue) lwidth(medium) ///
			|| scatter irr_beta_3 n if main==1, msymbol(O) mcolor(blue) ///
			title() legend(off) name(z_Mcaid_x_BigBoost90,replace)
		noisily disp "zMcaid_x_BigBoost"	
		gtoplevelsof sizemeasure if ub90<0
			
		*****************************************
		* Get the unique values of covariates
		*****************************************
		preserve
		capture drop fc2split 
		gen fc2split = stritrim(strtrim(fullcommand))
		split fc2split , parse(" ")
		drop fc2split
		reshape long fc2split, i(n) j(j)
			keep fc2split
			drop if inlist(fc2split," " , "," , ")" , "cluster(" , "cluster(CntyFE)" , "CntyFE" , "nbreg" , "covid_cases" , "z_Mcaid_x_BigBoost")
			bysort fc2split: keep if _n==1
			sort fc2split
			replace fc2split = subinstr(fc2split, ".", "_", .)
			levelsof fc2split ,clean
			global LEV `r(levels)'
		restore
	   foreach l in $LEV {
			noisily disp "`l'"
			  capture gen `l'  = strpos(fullcommand,"`l' ") >0 
	   }
	   foreach fe in CntyFE Beds20bins AvgResidents20bins {
			capture replace i_`fe' = 1 if  strpos(fullcommand,"`fe'")>0	
	   }
	   rename z_Beds Beds
	   rename z_Beds_2 Beds2
	   rename z_Beds_3 Beds3
	   rename i_Beds20bins Beds20bins
	   rename z_AvgResidents AvgResidents
	   rename z_AvgResidents_2 AvgResidents2
	   rename z_AvgResidents_3 AvgResidents3
	   rename i_AvgResidents20bins AvgResidents20bins
	   rename z_admissions_hundreds Admissions
	   rename Occupancy Occupancy
	   rename i_CntyFE OtherControls
		global V_LIST Beds Beds2 Beds3 Beds20bins AvgResidents AvgResidents2 AvgResidents3 AvgResidents20bins Admissions Occupancy OtherControls
		keep $V_LIST
	   *Save as matrix
	   mkmat  $V_LIST , matrix(VarPresent)
		   matrix Vt = VarPresent'
		   matrix list Vt
		heatplot Vt, colors(mono,reverse) legend(off) xlabel(,labsize(small)) ytitle("Included Variables",width(0)) ylabel(,  labsize(small)) ysize(2) xsize(6) name(vars, replace) plotregion(margin(zero)) graphregion(margin(t-5 b-12))

			graph combine z_Mcaid_x_BigBoost90 vars, cols(1) ysize(3) xsize(6) name(z_Mcaid_x_BigBoost,replace)
				graph drop z_Mcaid_x_BigBoost90 vars
